EOCSI EASI training session 3: Working with SAR data¶

This notebook will demonstrate how to use Sentinel-1 backscatter to detect the presence of water.

Introduction to Sentinel-1¶

Sentinel-1 carries a synthetic aperture radar (SAR) sensor, which can observe our planet even through cloud cover. Sentinel-1's SAR sensor is an active sensor; it emits light in the microwave range of the electormagnetic spectrum, and measures the amount that returns, known as backscatter. Smooth surfaces, like calm water, have very low backscatter - most of the light bounces away from the sensor (known as specular reflection). For rough surfaces, like dirt or vegetation, light is scattered in all directions, producing small backscatter signals (known as diffuse backscatter). Large, human-made structures, which feature both vertical and horizonal smooth surfaces, produce large backscatter signals (known as double bounce). As such, the intensity of Sentinel-1 backscatter can distinguish water from land, as shown in the image below:

image.png

In the image, the river appears dark (specular reflection), with the urban area appearing very bright (double bounce). For more information, see the article from GIS Geography on SAR.

This workshop¶

As of this workshop, Sentinel-1 backscatter is not yet available over South East Asia. Instead, this notebook will demonstrate how to load Sentinel-1 backscatter data over Africa, using Digital Earth Africa's STAC metadata. The Open Data Cube provides odc-stac, a useful library for loading data from STAC into an xarray, helping you work with it in the EASI platform. We will use this library in the workshop.

Important: If you want to work with Sentinel-1 as part of your project, please let us know and we can explore options to make data available to you.

The notebook contains the following steps:

  1. Import required packages for loading data and analysis
  2. Access Sentinel-1 analysis ready data from Digital Earth Africa
  3. Visualise data and export a multi-band, single time-slice GeoTIFF
  4. Applying speckle filter and converting the digital numbers to decibels (dB) values for analysis
  5. Use histogram analysis to determine the threshold for water classification
  6. Design a classifier to distinguish land and water
  7. Apply the classifier to the area of interest and interpret the results

Set up¶

Import required packages and functions¶

In [1]:
%matplotlib inline

# import required packages
import os

import matplotlib.pyplot as plt
import numpy as np
from dask.diagnostics import ProgressBar
from datacube.utils.cog import write_cog
from dea_tools.plotting import rgb, display_map
from odc.stac import configure_rio, stac_load
from pystac_client import Client
from scipy.ndimage import uniform_filter, variance
from skimage.filters import threshold_minimum

Make a directory in the tutorials folder to store results in¶

In [2]:
# make a results directory to store output

if not os.path.exists("results"):
    os.makedirs("results")

Configure s3 access¶

Digital Earth Africa data is stored on S3 in Cape Town, Africa. To load the data, we must configure rasterio with the appropriate AWS S3 endpoint. This is done with the odc.stac.configure_rio function.

The configuration below must be used when loading any Digital Earth Africa data through the STAC API.

In [3]:
# set AWS configuration with specific endpoint
configure_rio(
    cloud_defaults=True,
    aws={"aws_unsigned": True},
    AWS_S3_ENDPOINT="s3.af-south-1.amazonaws.com",
)

Open Digital Earth Africa's STAC endpoint¶

Be patient, as this will take a few minutes

In [4]:
# Open the stac catalogue
catalog = Client.open("https://explorer.digitalearth.africa/stac")
catalog
Out[4]:

Client: DEAfrica_data

ID: DEAfrica_data
Title: Digital Earth Africa
Description: Configure stac endpoint information in your Explorer `settings.env.py` file
type: Catalog
conformsTo: ['https://api.stacspec.org/v1.0.0-beta.2/core', 'https://api.stacspec.org/v1.0.0-beta.2/item-search']

Children

Only the first child shown

CollectionClient: Arrivals

ID: Arrivals
Title: Dataset Arrivals
Description: The most recently added Items to this index
title: Dataset Arrivals
type: Collection
properties: {}
providers: []
stac_extensions: []

Items

Only the first item shown

Item: 2e7917da-8db5-5380-84dc-5716ca447afa

ID: 2e7917da-8db5-5380-84dc-5716ca447afa
Bounding Box: [8.329440145578067, 37.842671646059465, 10.948502670082485, 39.95355032947056]
Datetime: 2022-09-29 10:00:59.906526+00:00
title: fc_ls_192033_2022-09-29
gsd: 30.0
created: 2022-10-05T03:48:52.037663Z
proj:epsg: 32632
platform: landsat-8
landsat:rmse: 7.097
odc:producer: digitalearthafrica.org
instruments: ['oli', 'tirs']
eo:cloud_cover: 39.58
view:sun_azimuth: 154.30324656
landsat:rmse_x: 4.804
landsat:rmse_y: 5.225
landsat:wrs_row: 33
odc:file_format: GeoTIFF
odc:region_code: 192033
constellation: Landsat
view:sun_elevation: 45.51499374
landsat:scene_id: LC81920332022272LGN00
landsat:wrs_path: 192
odc:product_family: fc
odc:dataset_version: 1.1.0
odc:naming_conventions: deafrica
landsat:collection_number: 2
landsat:collection_category: T1
proj:shape: [7821, 7691]
proj:transform: [30.0, 0.0, 439485.0, 0.0, -30.0, 4423215.0, 0.0, 0.0, 1.0]
datetime: 2022-09-29T10:00:59.906526Z
cubedash:region_code: 192033
stac_extensions: ['https://stac-extensions.github.io/eo/v1.0.0/schema.json', 'https://stac-extensions.github.io/projection/v1.0.0/schema.json', 'https://stac-extensions.github.io/view/v1.0.0/schema.json']

STAC Extensions

https://stac-extensions.github.io/eo/v1.0.0/schema.json
https://stac-extensions.github.io/projection/v1.0.0/schema.json
https://stac-extensions.github.io/view/v1.0.0/schema.json

Assets

Asset: bs

href: s3://deafrica-services/fc_ls/1-1-0/192/033/2022/09/29/fc_ls_192033_2022-09-29_bs.tif
Title: bs
Media type: image/tiff; application=geotiff; profile=cloud-optimized
Roles: ['data']
Owner:
eo:bands: [{'name': 'bs'}]
proj:epsg: 32632
proj:shape: [7821, 7691]
proj:transform: [30.0, 0.0, 439485.0, 0.0, -30.0, 4423215.0, 0.0, 0.0, 1.0]

Asset: pv

href: s3://deafrica-services/fc_ls/1-1-0/192/033/2022/09/29/fc_ls_192033_2022-09-29_pv.tif
Title: pv
Media type: image/tiff; application=geotiff; profile=cloud-optimized
Roles: ['data']
Owner:
eo:bands: [{'name': 'pv'}]
proj:epsg: 32632
proj:shape: [7821, 7691]
proj:transform: [30.0, 0.0, 439485.0, 0.0, -30.0, 4423215.0, 0.0, 0.0, 1.0]

Asset: ue

href: s3://deafrica-services/fc_ls/1-1-0/192/033/2022/09/29/fc_ls_192033_2022-09-29_ue.tif
Title: ue
Media type: image/tiff; application=geotiff; profile=cloud-optimized
Roles: ['data']
Owner:
eo:bands: [{'name': 'ue'}]
proj:epsg: 32632
proj:shape: [7821, 7691]
proj:transform: [30.0, 0.0, 439485.0, 0.0, -30.0, 4423215.0, 0.0, 0.0, 1.0]

Asset: npv

href: s3://deafrica-services/fc_ls/1-1-0/192/033/2022/09/29/fc_ls_192033_2022-09-29_npv.tif
Title: npv
Media type: image/tiff; application=geotiff; profile=cloud-optimized
Roles: ['data']
Owner:
eo:bands: [{'name': 'npv'}]
proj:epsg: 32632
proj:shape: [7821, 7691]
proj:transform: [30.0, 0.0, 439485.0, 0.0, -30.0, 4423215.0, 0.0, 0.0, 1.0]

Asset: Thumbnail image

href: s3://deafrica-services/fc_ls/1-1-0/192/033/2022/09/29/fc_ls_192033_2022-09-29_thumbnail.jpg
Title: Thumbnail image
Media type: image/jpeg
Roles: ['thumbnail']
Owner:

Asset: None

href: s3://deafrica-services/fc_ls/1-1-0/192/033/2022/09/29/fc_ls_192033_2022-09-29.sha1
Media type: text/plain
Roles: ['metadata']
Owner:

Asset: None

href: s3://deafrica-services/fc_ls/1-1-0/192/033/2022/09/29/fc_ls_192033_2022-09-29.proc-info.yaml
Media type: text/yaml
Roles: ['metadata']
Owner:

Links

Link:

Rel: self
Target: https://explorer.digitalearth.africa/stac/collections/fc_ls/items/2e7917da-8db5-5380-84dc-5716ca447afa
Media Type: application/json

Link:

ODC Dataset YAML

Rel: odc_yaml
Target: https://explorer.digitalearth.africa/dataset/2e7917da-8db5-5380-84dc-5716ca447afa.odc-metadata.yaml
Media Type: text/yaml

Link:

Rel: collection
Target: https://explorer.digitalearth.africa/stac/collections/fc_ls

Link:

ODC Product Overview

Rel: product_overview
Target: https://explorer.digitalearth.africa/product/fc_ls
Media Type: text/html

Link:

ODC Dataset Overview

Rel: alternative
Target: https://explorer.digitalearth.africa/dataset/2e7917da-8db5-5380-84dc-5716ca447afa
Media Type: text/html

Links

Link:

Rel: items
Target: https://explorer.digitalearth.africa/stac/arrivals/items

Link:

Digital Earth Africa

Rel: root
Target:
Media Type: application/json

Link:

Rel: self
Target: https://explorer.digitalearth.africa/stac/arrivals
Media Type: application/json

Link:

Digital Earth Africa

Rel: parent
Target:
Media Type: application/json

Items

Only the first item shown

Item: 2ba5be73-e8e7-5eed-84b9-b658e347e8e8

ID: 2ba5be73-e8e7-5eed-84b9-b658e347e8e8
Bounding Box: [-20.0, -40.000001192092896, 55.00000111758709, 40.0]
Datetime: 1981-01-01 00:00:00+00:00
title: chirps-v2.0_1981.01.01.stac-item
proj:epsg: 4326
proj:shape: [1600, 1500]
odc:product: rainfall_chirps_daily
proj:transform: [0.05000000074505806, 0.0, -20.0, 0.0, -0.05000000074505806, 40.0, 0.0, 0.0, 1.0]
odc:file_format: GeoTIFF
end_datetime: 1981-01-01T23:59:59Z
start_datetime: 1981-01-01T00:00:00Z
created: 2021-12-16T05:00:16.622284Z
datetime: 1981-01-01T00:00:00Z
cubedash:region_code: None
stac_extensions: ['https://stac-extensions.github.io/eo/v1.0.0/schema.json', 'https://stac-extensions.github.io/projection/v1.0.0/schema.json']

STAC Extensions

https://stac-extensions.github.io/eo/v1.0.0/schema.json
https://stac-extensions.github.io/projection/v1.0.0/schema.json

Assets

Asset: rainfall

href: s3://deafrica-input-datasets/rainfall_chirps_daily/1981/01/chirps-v2.0_1981.01.01.tif
Title: rainfall
Media type: image/tiff; application=geotiff; profile=cloud-optimized
Roles: ['data']
Owner:
eo:bands: [{'name': 'rainfall'}]
proj:epsg: 4326
proj:shape: [1600, 1500]
proj:transform: [0.05000000074505806, 0.0, -20.0, 0.0, -0.05000000074505806, 40.0, 0.0, 0.0, 1.0]

Links

Link:

Rel: self
Target: https://explorer.digitalearth.africa/stac/collections/rainfall_chirps_daily/items/2ba5be73-e8e7-5eed-84b9-b658e347e8e8
Media Type: application/json

Link:

ODC Dataset YAML

Rel: odc_yaml
Target: https://explorer.digitalearth.africa/dataset/2ba5be73-e8e7-5eed-84b9-b658e347e8e8.odc-metadata.yaml
Media Type: text/yaml

Link:

Rel: collection
Target: https://explorer.digitalearth.africa/stac/collections/rainfall_chirps_daily

Link:

ODC Product Overview

Rel: product_overview
Target: https://explorer.digitalearth.africa/product/rainfall_chirps_daily
Media Type: text/html

Link:

ODC Dataset Overview

Rel: alternative
Target: https://explorer.digitalearth.africa/dataset/2ba5be73-e8e7-5eed-84b9-b658e347e8e8
Media Type: text/html

Link:

Digital Earth Africa

Rel: root
Target:
Media Type: application/json

Links

Link:

Rel: self
Target: https://explorer.digitalearth.africa/stac
Media Type: application/json

Link:

Collections

Rel: children
Target: https://explorer.digitalearth.africa/stac/collections
Media Type: application/json
description: All product collections

Link:

Arrivals

Rel: child
Target:
Media Type: application/json
description: Most recently added items

Link:

Item Search

Rel: search
Target: https://explorer.digitalearth.africa/stac/search
Media Type: application/json

Link:

alos_palsar_mosaic

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/alos_palsar_mosaic
description: ALOS/PALSAR and ALOS-2/PALSAR-2 annual mosaic tiles generated for use in the Data Cube - 25m pixel spacing, WGS84. These tiles are derived from the orignal JAXA mosaics with conversion to GeoTIFF.

Link:

cci_landcover

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/cci_landcover
description: ESA Climate Change Initiative Land Cover

Link:

cgls_landcover

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/cgls_landcover
description: Copernicus Global Land Service, Land Use/Land Cover at 100 m

Link:

crop_mask

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/crop_mask
description: Annual cropland extent map produced by Digital Earth Africa.

Link:

crop_mask_central

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/crop_mask_central
description: Annual cropland extent map for Central Africa produced by Digital Earth Africa.

Link:

crop_mask_eastern

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/crop_mask_eastern
description: Annual cropland extent map for Eastern Africa produced by Digital Earth Africa.

Link:

crop_mask_indian_ocean

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/crop_mask_indian_ocean
description: Annual cropland extent map for Indian Ocean Africa produced by Digital Earth Africa.

Link:

crop_mask_northern

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/crop_mask_northern
description: Annual cropland extent map for Northern Africa produced by Digital Earth Africa.

Link:

crop_mask_sahel

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/crop_mask_sahel
description: Annual cropland extent map for Sahel Africa produced by Digital Earth Africa.

Link:

crop_mask_southeast

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/crop_mask_southeast
description: Annual cropland extent map for Southeast Africa produced by Digital Earth Africa.

Link:

crop_mask_southern

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/crop_mask_southern
description: Annual cropland extent map for Southern Africa produced by Digital Earth Africa.

Link:

crop_mask_western

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/crop_mask_western
description: Annual cropland extent map for Western Africa produced by Digital Earth Africa.

Link:

dem_cop_30

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/dem_cop_30
description: Copernicus DEM 30 m

Link:

dem_cop_90

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/dem_cop_90
description: Copernicus DEM 90 m

Link:

dem_srtm

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/dem_srtm
description: 1 second elevation model

Link:

dem_srtm_deriv

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/dem_srtm_deriv
description: 1 second elevation model derivatives

Link:

esa_worldcover

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/esa_worldcover
description: ESA World Cover, global 10 m land use/land cover data from 2020.

Link:

fc_ls

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/fc_ls
description: Landsat Fractional Cover Observations from Space

Link:

fc_ls_summary_annual

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/fc_ls_summary_annual
description: DE Africa Landsat Fractional Cover Percentiles

Link:

gm_ls5_ls7_annual

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/gm_ls5_ls7_annual
description: Surface Reflectance Annual Geometric Median and Median Absolute Deviations, Landsat 5 and Landsat 7

Link:

gm_ls5_ls7_annual_lowres

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/gm_ls5_ls7_annual_lowres
description: Surface Reflectance Annual Geometric Median and Median Absolute Deviations, Landsat 5 and Landsat 7. Low resolution version used for visualisations.

Link:

gm_ls8_annual

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/gm_ls8_annual
description: Surface Reflectance Annual Geometric Median and Median Absolute Deviations, Landsat 8

Link:

gm_ls8_annual_lowres

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/gm_ls8_annual_lowres
description: Surface Reflectance Annual Geometric Median and Median Absolute Deviations, Landsat 8. Low resolution version used for visualisations.

Link:

gm_ls8_ls9_annual

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/gm_ls8_ls9_annual
description: Surface Reflectance Annual Geometric Median and Median Absolute Deviations, Landsat 8 and Landsat 9

Link:

gm_ls8_ls9_annual_lowres

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/gm_ls8_ls9_annual_lowres
description: Surface Reflectance Annual Geometric Median and Median Absolute Deviations, Landsat 8 and Landsat 9. Low resolution version used for visualisations.

Link:

gm_s2_annual

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/gm_s2_annual
description: Surface Reflectance Annual Geometric Median and Median Absolute Deviations, Sentinel-2

Link:

gm_s2_annual_lowres

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/gm_s2_annual_lowres
description: Annual Geometric Median, Sentinel-2 - Low Resolution mosaic

Link:

gm_s2_semiannual

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/gm_s2_semiannual
description: Surface Reflectance Semiannual Geometric Median and Median Absolute Deviations, Sentinel-2

Link:

gm_s2_semiannual_lowres

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/gm_s2_semiannual_lowres
description: Surface Reflectance Semiannual Geometric Median and Median Absolute Deviations, Sentinel-2. Low resolution version used for visualisations.

Link:

gmw

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/gmw
description: Global Mangrove Watch data sourced from the UN Environment Program at https://data.unep-wcmc.org/datasets/45

Link:

io_lulc

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/io_lulc
description: Impact Observatory (ESRI) Landcover Classification 9 Class

Link:

jers_sar_mosaic

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/jers_sar_mosaic
description: JERS-1 SAR annual mosaic tiles generated for use in the Data Cube 25m pixel spacing, WGS84. These tiles are derived from the orignal JAXA mosaics with conversion to GeoTIFF.

Link:

ls5_sr

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/ls5_sr
description: USGS Landsat 5 Collection 2 Level-2 Surface Reflectance

Link:

ls5_st

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/ls5_st
description: USGS Landsat 5 Collection 2 Level-2 Surface Temperature

Link:

ls7_sr

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/ls7_sr
description: USGS Landsat 7 Collection 2 Level-2 Surface Reflectance

Link:

ls7_st

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/ls7_st
description: USGS Landsat 7 Collection 2 Level-2 Surface Temperature

Link:

ls8_sr

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/ls8_sr
description: USGS Landsat 8 Collection 2 Level-2 Surface Reflectance

Link:

ls8_st

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/ls8_st
description: USGS Landsat 8 Collection 2 Level-2 Surface Temperature

Link:

ls9_sr

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/ls9_sr
description: USGS Landsat 9 Collection 2 Level-2 Surface Reflectance

Link:

ls9_st

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/ls9_st
description: USGS Landsat 9 Collection 2 Level-2 Surface Temperature

Link:

nasadem

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/nasadem
description: NASADEM from Microsoft's Planetary Computer

Link:

ndvi_anomaly

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/ndvi_anomaly
description: Monthly NDVI Anomalies produced by Digital Earth Africa.

Link:

ndvi_climatology_ls

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/ndvi_climatology_ls
description: Monthly NDVI Climatologies produced by Digital Earth Africa.

Link:

pc_s2_annual

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/pc_s2_annual
description: Surface Reflectance Annual Clear Pixel Count, Sentinel-2

Link:

rainfall_chirps_daily

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/rainfall_chirps_daily
description: Rainfall Estimates from Rain Gauge and Satellite Observations

Link:

rainfall_chirps_monthly

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/rainfall_chirps_monthly
description: Rainfall Estimates from Rain Gauge and Satellite Observations

Link:

s1_rtc

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/s1_rtc
description: Sentinel 1 Gamma0 normalised radar backscatter

Link:

s2_l2a

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/s2_l2a
description: Sentinel-2a and Sentinel-2b imagery, processed to Level 2A (Surface Reflectance) and converted to Cloud Optimized GeoTIFFs

Link:

wofs_ls

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/wofs_ls
description: Historic Flood Mapping Water Observations from Space

Link:

wofs_ls_summary_alltime

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/wofs_ls_summary_alltime
description: Water Observations from Space Alltime Statistics

Link:

wofs_ls_summary_annual

Rel: child
Target: https://explorer.digitalearth.africa/stac/collections/wofs_ls_summary_annual
description: Water Observations from Space Annual Statistics

Link:

Digital Earth Africa

Rel: root
Target:
Media Type: application/json

Load Sentinel-1 from the Digital Earth Africa STAC endpoint¶

The specifics of Digital Earth Africa's Sentinel-1 backscatter product are available in the Digital Earth Africa Explorer

The data has five measurements, we will look at two today:

  • vv is the amount of vertically polarised backscatter received by the sensor
  • vh is the amount of horizontally polarised backscatter received by the sensor

Polarisation is the orientation of the light relative to the sensor, and different scattering surfaces can effect the backscatter polarisation. In the case of Sentinel-1, the sensor sends out vertically polarised light (v) and measures the amount of vertically polarised light (vv) or horizontally polarised light (vh) returned.

Steps for loading data¶

  1. Set the configuration
  2. Set the query
  3. Identify datasets that match the query
  4. Load the datasets

Set Collection Configuration¶

Digital Earth Africa's Sentinel-1 data (s1_rtc) needs to be configured with additional information about the data type and no-data values. This is provided in the dictionary below:

The configuration dictionary is determined from the product’s definition, available at https://explorer.digitalearth.africa/products/s1_rtc

In [5]:
# VV, VH, AREA = NaN, float32
# angle = 255, uint8
# mask = 0, uint8

# set the collection configuration, setting the products definition
config = {
    "s1_rtc": {
        "assets": {
            "*": {
                "data_type": "float32",
                "nodata": float("nan"),
                "unit": "1",
            },
            "angle": {
                "data_type": "uint8",
                "nodata": 255,
                "unit": "1",
            },
            "mask": {
                "data_type": "uint8",
                "nodata": 0,
                "unit": "1",
            },
        },
    }
}

Define Query Parameters¶

For this example, we'll use a bounding box, but you could come back an update this to use a polygon by following training session 2

When working with STAC, we use the word "collection" to refer to a dataset, rather than "product" when working with the ODC. For more information on the differences, see the odc-stac documentation

In [6]:
# Set a bounding box using latitude and longitude coordinates
# [xmin, ymin, xmax, ymax]
bbox = [-16.34, 12.5699, -16.24, 12.67]

# timeframe
single_date = "2020"

# Set the STAC collections
collections = ["s1_rtc"]
In [7]:
# View area of interest
display_map(x=(bbox[0], bbox[2]), y=(bbox[1],bbox[3]))
Out[7]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Identify datasets that match the query¶

The connection to Digital Earth Africa's STAC endpoint provides a search function which returns dataset ids matching the parameters of the search. Data can then be loaded by accessing the corresponding files for each id.

In [8]:
# Build a query with the set parameters
query = catalog.search(bbox=bbox, collections=collections, datetime=single_date)

# Search the STAC catalog for all items matching the query
items = list(query.get_items())
print(f"Found: {len(items):d} datasets")
Found: 31 datasets

Load the Data¶

In this section we will use the stac_load function to only load data within the defined parameters set above.

This data will be lazy-loaded with dask, meaning the it will not be loaded into memory until required, such as when it is displayed.

In [9]:
crs = "EPSG:6933"
resolution = 20

ds = stac_load(
    items,
    bands=[
        "vv",
        "vh",
    ],
    crs=crs,
    resolution=resolution,
    chunks={},
    groupby="solar_day",
    stac_cfg=config,
    bbox=bbox,
)

# View the Xarray Dataset - loads as dask array
ds
Out[9]:
<xarray.Dataset>
Dimensions:      (y: 624, x: 484, time: 31)
Coordinates:
  * y            (y) float64 1.604e+06 1.604e+06 ... 1.591e+06 1.591e+06
  * x            (x) float64 -1.577e+06 -1.577e+06 ... -1.567e+06 -1.567e+06
    spatial_ref  int32 6933
  * time         (time) datetime64[ns] 2020-01-04T19:17:34.083167 ... 2020-12...
Data variables:
    vv           (time, y, x) float32 dask.array<chunksize=(1, 624, 484), meta=np.ndarray>
    vh           (time, y, x) float32 dask.array<chunksize=(1, 624, 484), meta=np.ndarray>
xarray.Dataset
    • y: 624
    • x: 484
    • time: 31
    • y
      (y)
      float64
      1.604e+06 1.604e+06 ... 1.591e+06
      units :
      metre
      resolution :
      -20.0
      crs :
      EPSG:6933
      array([1603550., 1603530., 1603510., ..., 1591130., 1591110., 1591090.])
    • x
      (x)
      float64
      -1.577e+06 ... -1.567e+06
      units :
      metre
      resolution :
      20.0
      crs :
      EPSG:6933
      array([-1576590., -1576570., -1576550., ..., -1566970., -1566950., -1566930.])
    • spatial_ref
      ()
      int32
      6933
      spatial_ref :
      PROJCRS["WGS 84 / NSIDC EASE-Grid 2.0 Global",BASEGEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]],CONVERSION["US NSIDC EASE-Grid 2.0 Global",METHOD["Lambert Cylindrical Equal Area",ID["EPSG",9835]],PARAMETER["Latitude of 1st standard parallel",30,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting (X)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["northing (Y)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Environmental science - used as basis for EASE grid."],AREA["World between 86°S and 86°N."],BBOX[-86,-180,86,180]],ID["EPSG",6933]]
      crs_wkt :
      PROJCRS["WGS 84 / NSIDC EASE-Grid 2.0 Global",BASEGEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]],CONVERSION["US NSIDC EASE-Grid 2.0 Global",METHOD["Lambert Cylindrical Equal Area",ID["EPSG",9835]],PARAMETER["Latitude of 1st standard parallel",30,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting (X)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["northing (Y)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Environmental science - used as basis for EASE grid."],AREA["World between 86°S and 86°N."],BBOX[-86,-180,86,180]],ID["EPSG",6933]]
      semi_major_axis :
      6378137.0
      semi_minor_axis :
      6356752.314245179
      inverse_flattening :
      298.257223563
      reference_ellipsoid_name :
      WGS 84
      longitude_of_prime_meridian :
      0.0
      prime_meridian_name :
      Greenwich
      geographic_crs_name :
      WGS 84
      horizontal_datum_name :
      World Geodetic System 1984 ensemble
      projected_crs_name :
      WGS 84 / NSIDC EASE-Grid 2.0 Global
      grid_mapping_name :
      lambert_cylindrical_equal_area
      standard_parallel :
      30.0
      longitude_of_central_meridian :
      0.0
      false_easting :
      0.0
      false_northing :
      0.0
      array(6933, dtype=int32)
    • time
      (time)
      datetime64[ns]
      2020-01-04T19:17:34.083167 ... 2...
      array(['2020-01-04T19:17:34.083167000', '2020-01-16T19:17:33.556635000',
             '2020-01-28T19:17:33.185197000', '2020-02-09T19:17:32.895481000',
             '2020-02-21T19:17:32.683429000', '2020-03-04T19:17:32.717588000',
             '2020-03-16T19:17:32.881860000', '2020-03-28T19:17:32.985420000',
             '2020-04-09T19:17:33.213406000', '2020-04-21T19:17:33.757711000',
             '2020-05-03T19:17:34.346095000', '2020-05-15T19:17:35.068275000',
             '2020-05-27T19:17:35.817349000', '2020-06-08T19:17:36.497360000',
             '2020-06-20T19:17:37.233415000', '2020-07-02T19:17:37.706227000',
             '2020-07-14T19:17:38.462335000', '2020-07-26T19:17:39.446203000',
             '2020-08-07T19:17:39.987574000', '2020-08-19T19:17:40.585644000',
             '2020-08-31T19:17:41.407335000', '2020-09-12T19:17:41.863551000',
             '2020-09-24T19:17:42.242509000', '2020-10-06T19:17:42.525580000',
             '2020-10-18T19:17:42.548824000', '2020-10-30T19:17:42.569753000',
             '2020-11-11T19:17:42.271800000', '2020-11-23T19:17:41.906190000',
             '2020-12-05T19:17:41.558699000', '2020-12-17T19:17:41.078736000',
             '2020-12-29T19:17:40.517695000'], dtype='datetime64[ns]')
    • vv
      (time, y, x)
      float32
      dask.array<chunksize=(1, 624, 484), meta=np.ndarray>
      nodata :
      nan
      Array Chunk
      Bytes 35.72 MiB 1.15 MiB
      Shape (31, 624, 484) (1, 624, 484)
      Count 1 Graph Layer 31 Chunks
      Type float32 numpy.ndarray
      484 624 31
    • vh
      (time, y, x)
      float32
      dask.array<chunksize=(1, 624, 484), meta=np.ndarray>
      nodata :
      nan
      Array Chunk
      Bytes 35.72 MiB 1.15 MiB
      Shape (31, 624, 484) (1, 624, 484)
      Count 1 Graph Layer 31 Chunks
      Type float32 numpy.ndarray
      484 624 31
In [10]:
# Load the data in memory using dask's .compute() function

with ProgressBar():
    ds = ds.compute()
[########################################] | 100% Completed | 11.15 s

Visualise the time series¶

The vv and vh bands can be viewed individually, or as a RGB style image through combination with their ratio as the third band.

In [11]:
# Selecting a few images from the loaded S1 to visualise
timesteps = [2, 4, 6, 9, 11]
In [12]:
# Plot VV polarisation for specific timeframe
ds.vv.isel(time=timesteps).plot(cmap="Greys_r", robust=True, col="time", col_wrap=5);
In [13]:
# Plot VH polarisation for specific timeframe
ds.vh.isel(time=timesteps).plot(cmap="Greys_r", robust=True, col="time", col_wrap=5);

For the RGB visualization below, the ratio between VH and VV is added as a third measurement band.

In [14]:
# VH/VV is a potentially useful third feature after VV and VH
ds["vh/vv"] = ds.vh / ds.vv
In [15]:
# median values are used to scale the measurements so they have a similar range for visualization
med_s1 = ds.median()
In [16]:
# plotting an RGB image for selected timesteps
rgb(
    ds[["vv", "vh", "vh/vv"]] / med_s1,
    bands=["vv", "vh", "vh/vv"],
    index=timesteps,
    col_wrap=5,
);

Export a multi-band, single time-slice GeoTIFF¶

The next three cells will export the first observation (time=0) to a GeoTiff.

If you want to export multiple observations, refer to training session 2.

In [17]:
da = ds.to_array()
In [18]:
rgb_ga = da.isel(time=0)
In [19]:
write_cog(geo_im=rgb_ga, fname="results/sar_rgb.tif", overwrite=True)
Out[19]:
PosixPath('results/sar_rgb.tif')

Apply speckle filtering¶

A common processing step when working with SAR data is to remove , which causes SAR data to have a grainy appearence. Radar observations appear speckly) due to random interference of coherent signals from target scatters. The speckle noise can be reduced by averaging pixel values over an area or over time. However, averaging over a fixed window smoothes out real local spatial variation and leads to reduced spatial resolution. An adaptive approach that takes into account local homogeneity is therefore preferred.

Below, we display the unfiltered data, then apply the Lee filter, one of the popular adaptive speckle filters.

In [20]:
# images appear smoother after speckle filtering
fig, ax = plt.subplots(1, 2, figsize=(15, 5))
ds["vv"].isel(time=3).plot(ax=ax[0], robust=True)
ds["vh"].isel(time=3).plot(ax=ax[1], robust=True)
ax[0].set_title("VV")
ax[1].set_title("VH")
plt.tight_layout();
In [21]:
# defining a function to apply lee filtering on S1 image
def lee_filter(da, size):
    """
    Apply lee filter of specified window size.
    Adapted from https://stackoverflow.com/questions/39785970/speckle-lee-filter-in-python

    """
    img = da.values
    img_mean = uniform_filter(img, (size, size))
    img_sqr_mean = uniform_filter(img**2, (size, size))
    img_variance = img_sqr_mean - img_mean**2

    overall_variance = variance(img)

    img_weights = img_variance / (img_variance + overall_variance)
    img_output = img_mean + img_weights * (img - img_mean)

    return img_output

Now that we’ve defined the filter, we can run it on the VV and VH data. You might have noticed that the function takes a size argument. This will change how blurred the image becomes after smoothing. We’ve picked a default value for this analysis, but you can experiement with this if you’re interested.

In [22]:
# The lee filter above doesn't handle null values
# We therefore set null values to 0 before applying the filter
valid = ds.where(~np.isinf(ds))
ds = ds.where(valid, 0)

# Create a new entry in dataset corresponding to filtered VV and VH data
ds["filtered_vv"] = ds.vv.groupby("time").apply(lee_filter, size=7)
ds["filtered_vh"] = ds.vh.groupby("time").apply(lee_filter, size=7)

# Null pixels should remain null
ds["filtered_vv"] = ds.filtered_vv.where(valid.vv, np.nan)
ds["filtered_vh"] = ds.filtered_vh.where(valid.vh, np.nan)
In [23]:
# images appear smoother after speckle filtering
fig, ax = plt.subplots(1, 2, figsize=(15, 5))
ds["filtered_vv"].isel(time=3).plot(ax=ax[0], robust=True)
ds["filtered_vh"].isel(time=3).plot(ax=ax[1], robust=True)
ax[0].set_title("VV")
ax[1].set_title("VH")
plt.tight_layout();

Designing a threshold-based water classifier¶

Similar to training session 2, it is possible to view the histograms for the vv and vh bands and identify an automatic threshold that seperates two key features in the landscape: water, which has very low backscatter, and land, which has higher backscatter.

When working with radar backscatter, it is common to work with the data in units of decibels (dB), rather than digital number (DN) as measured by the satellite. To convert from DN to dB, we use the following formula:

$$\text{dB} = 10 \times \log_{10}(\text{DN})$$
In [24]:
# convert digital numbers to dB
ds["filtered_vv"] = 10 * np.log10(ds.filtered_vv)
ds["filtered_vh"] = 10 * np.log10(ds.filtered_vh)
In [25]:
# histogram analysis for S1
fig = plt.figure(figsize=(12, 3))
ds.filtered_vh.plot.hist(bins=1000, label="VH filtered")
ds.filtered_vv.plot.hist(bins=1000, label="VV filtered", alpha=0.5)
plt.xlim(-40, -1)
plt.legend()
plt.xlabel("DN values in(dB)")
plt.title("Comparison of Lee filtered VH and VV polarisation bands");

The vv band has a clear separation between the two classes. There is a smaller number of low value pixels, peaking around -30 dB (likely water), and a larger number of high value pixels, peaking around -15 dB (likely land).

Build and apply the water classifier¶

The histogram for VH backscatter shows a bimodal distribution with low values over water and high values over land. The VV histogram has multiple peaks and less obvious seperation between water and land.

We therefore build a classifier based on VH backscatter. We choose a threshold to separate land and water: pixels with values below the threshold are water, and pixels with values above the threshold are not water (land).

There are several ways to determine the threshold. Here, we use the threshod_minimum function implemented in the skimage package to determine the threshold from the VH histogram automatically. This method computes the histogram for all backscatter values, smooths it until there are only two maxima and find the minimum in between as the threshold.

In [26]:
# Convert to numpy ndarray, remove nan values, and calculate automatic threshold
vh_new = ds["filtered_vh"].values
vh_new = vh_new[~np.isnan(vh_new)]

threshold_vh = threshold_minimum(vh_new)

print(threshold_vh)
-25.661585
In [27]:
# visualise the threshold
fig, ax = plt.subplots(figsize=(15, 3))
ds.filtered_vh.plot.hist(bins=1000, label="VH filtered")
plt.xlim(-40, -5)
ax.axvspan(xmin=-40.0, xmax=threshold_vh, alpha=0.25, color="green", label="Water")
ax.axvspan(xmin=threshold_vh, xmax=-5, alpha=0.25, color="red", label="Not Water")
plt.legend()
plt.xlabel("VH (dB)")
plt.title("Effect of the classifier")
plt.show()
In [28]:
# define the classifier, values lower than the threshold are water pixels, return dataset containing water pixels
def S1_water_classifier(da, threshold=threshold_vh):
    water_data_array = da < threshold
    return water_data_array.to_dataset(name="s1_water")
In [29]:
# return classified data product
ds["water"] = S1_water_classifier(ds.filtered_vh).s1_water
ds
Out[29]:
<xarray.Dataset>
Dimensions:      (time: 31, y: 624, x: 484)
Coordinates:
  * y            (y) float64 1.604e+06 1.604e+06 ... 1.591e+06 1.591e+06
  * x            (x) float64 -1.577e+06 -1.577e+06 ... -1.567e+06 -1.567e+06
    spatial_ref  int32 6933
  * time         (time) datetime64[ns] 2020-01-04T19:17:34.083167 ... 2020-12...
Data variables:
    vv           (time, y, x) float32 0.05976 0.04359 0.02704 ... 0.2923 0.1862
    vh           (time, y, x) float32 0.005369 0.004237 ... 0.01821 0.01918
    vh/vv        (time, y, x) float32 0.08984 0.09722 0.1888 ... 0.06229 0.103
    filtered_vv  (time, y, x) float32 -12.21 -12.32 -12.53 ... -7.404 -7.758
    filtered_vh  (time, y, x) float32 -22.11 -22.5 -22.74 ... -19.43 -19.19
    water        (time, y, x) bool False False False False ... False False False
xarray.Dataset
    • time: 31
    • y: 624
    • x: 484
    • y
      (y)
      float64
      1.604e+06 1.604e+06 ... 1.591e+06
      units :
      metre
      resolution :
      -20.0
      crs :
      EPSG:6933
      array([1603550., 1603530., 1603510., ..., 1591130., 1591110., 1591090.])
    • x
      (x)
      float64
      -1.577e+06 ... -1.567e+06
      units :
      metre
      resolution :
      20.0
      crs :
      EPSG:6933
      array([-1576590., -1576570., -1576550., ..., -1566970., -1566950., -1566930.])
    • spatial_ref
      ()
      int32
      6933
      spatial_ref :
      PROJCRS["WGS 84 / NSIDC EASE-Grid 2.0 Global",BASEGEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]],CONVERSION["US NSIDC EASE-Grid 2.0 Global",METHOD["Lambert Cylindrical Equal Area",ID["EPSG",9835]],PARAMETER["Latitude of 1st standard parallel",30,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting (X)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["northing (Y)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Environmental science - used as basis for EASE grid."],AREA["World between 86°S and 86°N."],BBOX[-86,-180,86,180]],ID["EPSG",6933]]
      crs_wkt :
      PROJCRS["WGS 84 / NSIDC EASE-Grid 2.0 Global",BASEGEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]],CONVERSION["US NSIDC EASE-Grid 2.0 Global",METHOD["Lambert Cylindrical Equal Area",ID["EPSG",9835]],PARAMETER["Latitude of 1st standard parallel",30,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting (X)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["northing (Y)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Environmental science - used as basis for EASE grid."],AREA["World between 86°S and 86°N."],BBOX[-86,-180,86,180]],ID["EPSG",6933]]
      semi_major_axis :
      6378137.0
      semi_minor_axis :
      6356752.314245179
      inverse_flattening :
      298.257223563
      reference_ellipsoid_name :
      WGS 84
      longitude_of_prime_meridian :
      0.0
      prime_meridian_name :
      Greenwich
      geographic_crs_name :
      WGS 84
      horizontal_datum_name :
      World Geodetic System 1984 ensemble
      projected_crs_name :
      WGS 84 / NSIDC EASE-Grid 2.0 Global
      grid_mapping_name :
      lambert_cylindrical_equal_area
      standard_parallel :
      30.0
      longitude_of_central_meridian :
      0.0
      false_easting :
      0.0
      false_northing :
      0.0
      array(6933, dtype=int32)
    • time
      (time)
      datetime64[ns]
      2020-01-04T19:17:34.083167 ... 2...
      array(['2020-01-04T19:17:34.083167000', '2020-01-16T19:17:33.556635000',
             '2020-01-28T19:17:33.185197000', '2020-02-09T19:17:32.895481000',
             '2020-02-21T19:17:32.683429000', '2020-03-04T19:17:32.717588000',
             '2020-03-16T19:17:32.881860000', '2020-03-28T19:17:32.985420000',
             '2020-04-09T19:17:33.213406000', '2020-04-21T19:17:33.757711000',
             '2020-05-03T19:17:34.346095000', '2020-05-15T19:17:35.068275000',
             '2020-05-27T19:17:35.817349000', '2020-06-08T19:17:36.497360000',
             '2020-06-20T19:17:37.233415000', '2020-07-02T19:17:37.706227000',
             '2020-07-14T19:17:38.462335000', '2020-07-26T19:17:39.446203000',
             '2020-08-07T19:17:39.987574000', '2020-08-19T19:17:40.585644000',
             '2020-08-31T19:17:41.407335000', '2020-09-12T19:17:41.863551000',
             '2020-09-24T19:17:42.242509000', '2020-10-06T19:17:42.525580000',
             '2020-10-18T19:17:42.548824000', '2020-10-30T19:17:42.569753000',
             '2020-11-11T19:17:42.271800000', '2020-11-23T19:17:41.906190000',
             '2020-12-05T19:17:41.558699000', '2020-12-17T19:17:41.078736000',
             '2020-12-29T19:17:40.517695000'], dtype='datetime64[ns]')
    • vv
      (time, y, x)
      float32
      0.05976 0.04359 ... 0.2923 0.1862
      nodata :
      nan
      array([[[0.05975673, 0.0435871 , 0.0270416 , ..., 0.14097616,
               0.05979058, 0.0423508 ],
              [0.07631498, 0.07984253, 0.04585955, ..., 0.10154395,
               0.08561386, 0.05674401],
              [0.07631498, 0.07984253, 0.04585955, ..., 0.10154395,
               0.08561386, 0.05674401],
              ...,
              [0.0089145 , 0.00657869, 0.00835743, ..., 0.07020288,
               0.0757059 , 0.07646044],
              [0.00804577, 0.00638901, 0.00634187, ..., 0.0840879 ,
               0.11501899, 0.0833811 ],
              [0.00823483, 0.00586537, 0.00580274, ..., 0.13835531,
               0.13704413, 0.08174172]],
      
             [[0.05852966, 0.03650849, 0.02030807, ..., 0.1308103 ,
               0.05439572, 0.09825845],
              [0.07354559, 0.0284885 , 0.02448474, ..., 0.08454306,
               0.08480984, 0.10215209],
              [0.07354559, 0.0284885 , 0.02448474, ..., 0.08454306,
               0.08480984, 0.10215209],
      ...
              [0.06990052, 0.05567687, 0.06165422, ..., 0.04058722,
               0.07384066, 0.07156248],
              [0.10927196, 0.11126049, 0.08923999, ..., 0.0606725 ,
               0.09735471, 0.09874935],
              [0.19658579, 0.14020544, 0.08757   , ..., 0.09544989,
               0.13431168, 0.10657003]],
      
             [[0.06696489, 0.05246656, 0.03226287, ..., 0.0760254 ,
               0.08146016, 0.08312523],
              [0.07243357, 0.06536565, 0.04068712, ..., 0.06209026,
               0.08140337, 0.07891082],
              [0.07243357, 0.06536565, 0.04068712, ..., 0.06209026,
               0.08140337, 0.07891082],
              ...,
              [0.12821572, 0.08549089, 0.07939965, ..., 0.0908902 ,
               0.16557105, 0.16005161],
              [0.1490065 , 0.09286911, 0.07969621, ..., 0.1496388 ,
               0.18705945, 0.15205307],
              [0.12625034, 0.13715085, 0.12770365, ..., 0.22266619,
               0.29232228, 0.18619347]]], dtype=float32)
    • vh
      (time, y, x)
      float32
      0.005369 0.004237 ... 0.01918
      nodata :
      nan
      array([[[0.00536875, 0.00423742, 0.00510582, ..., 0.01670616,
               0.01049374, 0.01557176],
              [0.0076937 , 0.00664391, 0.00559766, ..., 0.0151424 ,
               0.00927143, 0.00849948],
              [0.0076937 , 0.00664391, 0.00559766, ..., 0.0151424 ,
               0.00927143, 0.00849948],
              ...,
              [0.00100636, 0.00095855, 0.00077952, ..., 0.00790127,
               0.00672355, 0.00669588],
              [0.00134101, 0.00109705, 0.00218619, ..., 0.00701612,
               0.00887922, 0.01061779],
              [0.00438691, 0.0015797 , 0.00294775, ..., 0.00504207,
               0.00683029, 0.01434511]],
      
             [[0.00805709, 0.0017193 , 0.00372427, ..., 0.01851185,
               0.01328646, 0.02111009],
              [0.0171579 , 0.0074999 , 0.00185625, ..., 0.02620168,
               0.01923757, 0.02179751],
              [0.0171579 , 0.0074999 , 0.00185625, ..., 0.02620168,
               0.01923757, 0.02179751],
      ...
              [0.00309067, 0.00181976, 0.00131528, ..., 0.00295775,
               0.00671206, 0.01157954],
              [0.00639177, 0.0037103 , 0.00299319, ..., 0.00447364,
               0.01041666, 0.01462556],
              [0.00500921, 0.00333671, 0.0055453 , ..., 0.00662705,
               0.00965548, 0.01567301]],
      
             [[0.01181371, 0.00454491, 0.00296509, ..., 0.02196231,
               0.02373657, 0.03063028],
              [0.01960026, 0.00905047, 0.00537738, ..., 0.01216274,
               0.01950232, 0.01728071],
              [0.01960026, 0.00905047, 0.00537738, ..., 0.01216274,
               0.01950232, 0.01728071],
              ...,
              [0.0027416 , 0.0014938 , 0.00207409, ..., 0.00686023,
               0.01253865, 0.0127871 ],
              [0.00795899, 0.00526682, 0.00277987, ..., 0.00700331,
               0.0157548 , 0.01733486],
              [0.00934836, 0.00607872, 0.00435709, ..., 0.00987481,
               0.01820743, 0.01917654]]], dtype=float32)
    • vh/vv
      (time, y, x)
      float32
      0.08984 0.09722 ... 0.06229 0.103
      array([[[0.08984352, 0.09721721, 0.18881354, ..., 0.11850342,
               0.1755083 , 0.36768514],
              [0.1008151 , 0.08321262, 0.12206098, ..., 0.14912167,
               0.10829357, 0.14978632],
              [0.1008151 , 0.08321262, 0.12206098, ..., 0.14912167,
               0.10829357, 0.14978632],
              ...,
              [0.11289068, 0.14570576, 0.09327213, ..., 0.11254911,
               0.0888115 , 0.08757319],
              [0.16667314, 0.1717087 , 0.34472248, ..., 0.08343796,
               0.07719785, 0.12734051],
              [0.53272575, 0.26932558, 0.507993  , ..., 0.03644292,
               0.04984008, 0.17549309]],
      
             [[0.1376582 , 0.04709305, 0.18338849, ..., 0.14151673,
               0.24425565, 0.21484244],
              [0.2332961 , 0.26326066, 0.0758126 , ..., 0.3099211 ,
               0.22683184, 0.21338288],
              [0.2332961 , 0.26326066, 0.0758126 , ..., 0.3099211 ,
               0.22683184, 0.21338288],
      ...
              [0.04421525, 0.03268425, 0.02133309, ..., 0.07287388,
               0.0908992 , 0.16181019],
              [0.05849414, 0.03334785, 0.0335409 , ..., 0.07373427,
               0.10699694, 0.1481079 ],
              [0.02548104, 0.02379869, 0.06332423, ..., 0.06942965,
               0.0718886 , 0.14706774]],
      
             [[0.17641653, 0.0866248 , 0.09190405, ..., 0.28888127,
               0.29138872, 0.36848354],
              [0.2705964 , 0.13845904, 0.1321642 , ..., 0.19588806,
               0.23957634, 0.21899039],
              [0.2705964 , 0.13845904, 0.1321642 , ..., 0.19588806,
               0.23957634, 0.21899039],
              ...,
              [0.02138271, 0.01747317, 0.02612215, ..., 0.07547821,
               0.07572971, 0.0798936 ],
              [0.05341373, 0.05671225, 0.03488087, ..., 0.0468014 ,
               0.08422349, 0.11400531],
              [0.07404618, 0.04432142, 0.03411876, ..., 0.04434805,
               0.06228545, 0.10299257]]], dtype=float32)
    • filtered_vv
      (time, y, x)
      float32
      -12.21 -12.32 ... -7.404 -7.758
      array([[[-12.209671 , -12.318764 , -12.5348015, ...,  -9.409913 ,
               -10.068176 , -10.280441 ],
              [-11.947147 , -12.076376 , -12.435654 , ...,  -9.80544  ,
               -10.10659  , -10.307478 ],
              [-11.758413 , -11.978493 , -12.411015 , ..., -10.039537 ,
               -10.319677 , -10.482844 ],
              ...,
              [-16.77635  , -15.947324 , -14.539786 , ..., -10.342581 ,
               -10.526849 , -10.513917 ],
              [-18.261837 , -17.146393 , -15.378963 , ..., -10.233499 ,
               -10.401096 , -10.471428 ],
              [-19.406406 , -18.0442   , -15.906037 , ..., -10.129518 ,
               -10.323764 , -10.425333 ]],
      
             [[-13.622372 , -13.469177 , -13.384103 , ...,  -8.702457 ,
                -9.500488 ,  -9.745563 ],
              [-13.313275 , -13.2765665, -13.2574835, ...,  -9.080205 ,
                -9.4850025,  -9.734879 ],
              [-13.005913 , -13.0680485, -13.180083 , ...,  -9.412878 ,
                -9.797277 ,  -9.978608 ],
      ...
              [-10.222945 , -10.136291 ,  -9.742037 , ..., -11.069022 ,
               -10.593851 , -10.433159 ],
              [-10.044184 ,  -9.753965 ,  -9.492738 , ..., -11.156181 ,
               -10.721802 , -10.608507 ],
              [ -9.926697 ,  -9.670756 ,  -9.494183 , ..., -11.239136 ,
               -10.844025 , -10.814868 ]],
      
             [[-12.57022  , -12.558077 , -12.593336 , ...,  -9.435576 ,
                -9.828341 , -10.46761  ],
              [-12.542982 , -12.543735 , -12.645113 , ...,  -9.712024 ,
                -9.999484 , -10.620775 ],
              [-12.584721 , -12.63774  , -12.782244 , ...,  -9.964315 ,
               -10.255446 , -10.811843 ],
              ...,
              [ -9.5582695,  -9.512032 ,  -9.325964 , ...,  -9.083597 ,
                -8.057753 ,  -7.933816 ],
              [ -9.448257 ,  -9.346028 ,  -9.140655 , ...,  -8.610491 ,
                -7.8879175,  -7.8889823],
              [ -9.492592 ,  -9.271316 ,  -8.9779005, ...,  -8.2096615,
                -7.403696 ,  -7.7577324]]], dtype=float32)
    • filtered_vh
      (time, y, x)
      float32
      -22.11 -22.5 ... -19.43 -19.19
      array([[[-22.10778 , -22.499306, -22.742874, ..., -16.517765,
               -17.630348, -18.151491],
              [-21.425676, -21.754707, -22.119017, ..., -16.674692,
               -17.658905, -18.23654 ],
              [-21.055046, -21.393196, -21.769602, ..., -16.873058,
               -17.757967, -18.242924],
              ...,
              [-26.909714, -25.510849, -23.160189, ..., -21.000069,
               -21.096561, -21.103098],
              [-27.298307, -25.668163, -22.916256, ..., -20.909325,
               -21.037485, -21.121042],
              [-27.485008, -25.724113, -22.893497, ..., -20.864985,
               -21.00068 , -21.037996]],
      
             [[-21.289045, -21.482616, -21.281376, ..., -16.34781 ,
               -16.690834, -16.682783],
              [-20.906261, -21.232002, -21.350374, ..., -16.323765,
               -16.779701, -16.854565],
              [-20.460081, -20.938042, -21.165398, ..., -16.554628,
               -17.074175, -17.178465],
      ...
              [-24.387226, -23.604277, -22.222422, ..., -21.363308,
               -21.253132, -21.20169 ],
              [-24.129265, -23.144024, -21.7937  , ..., -21.018534,
               -20.832333, -20.771214],
              [-24.35956 , -23.076763, -21.523916, ..., -21.005556,
               -20.863659, -20.758486]],
      
             [[-20.333673, -20.740917, -20.705425, ..., -15.743235,
               -15.869292, -15.967106],
              [-20.385427, -20.811905, -20.875557, ..., -16.370245,
               -16.21775 , -16.752844],
              [-20.655857, -21.03719 , -21.1148  , ..., -16.605793,
               -16.398638, -16.86529 ],
              ...,
              [-23.113844, -22.357132, -20.71432 , ..., -20.521206,
               -19.71997 , -19.462969],
              [-23.116661, -22.170828, -20.456417, ..., -20.318947,
               -19.438139, -19.181206],
              [-23.38486 , -22.341301, -20.450977, ..., -20.277924,
               -19.427477, -19.187847]]], dtype=float32)
    • water
      (time, y, x)
      bool
      False False False ... False False
      array([[[False, False, False, ..., False, False, False],
              [False, False, False, ..., False, False, False],
              [False, False, False, ..., False, False, False],
              ...,
              [ True, False, False, ..., False, False, False],
              [ True,  True, False, ..., False, False, False],
              [ True,  True, False, ..., False, False, False]],
      
             [[False, False, False, ..., False, False, False],
              [False, False, False, ..., False, False, False],
              [False, False, False, ..., False, False, False],
              ...,
              [ True,  True, False, ..., False, False, False],
              [ True,  True, False, ..., False, False, False],
              [ True,  True, False, ..., False, False, False]],
      
             [[False, False, False, ..., False, False, False],
              [False, False, False, ..., False, False, False],
              [False, False, False, ..., False, False, False],
              ...,
      ...
              ...,
              [False, False, False, ..., False, False, False],
              [ True, False, False, ..., False, False, False],
              [ True, False, False, ..., False, False, False]],
      
             [[False, False, False, ..., False, False, False],
              [False, False, False, ..., False, False, False],
              [False, False, False, ..., False, False, False],
              ...,
              [False, False, False, ..., False, False, False],
              [False, False, False, ..., False, False, False],
              [False, False, False, ..., False, False, False]],
      
             [[False, False, False, ..., False, False, False],
              [False, False, False, ..., False, False, False],
              [False, False, False, ..., False, False, False],
              ...,
              [False, False, False, ..., False, False, False],
              [False, False, False, ..., False, False, False],
              [False, False, False, ..., False, False, False]]])

Assessment with mean¶

We can now view the image with our classification. The classifier returns either True or False for each pixel. To investigate the landscape, we want to check which pixels are always water and which are always land. Conveniently, Python encodes True = 1 and False = 0. If we plot the average classified pixel value, pixels that are always water will have an average value of 1 and pixels that are always land will have an average of 0. Pixels that are sometimes water and sometimes land will have an average between these values.

The following cell plots the average classified pixel value over time. What can you learn about the landscape from this summary?

In [30]:
# Plot the mean of each classified pixel value
plt.figure(figsize=(8, 7))
ds.water.mean(dim="time").plot(cmap="RdBu")
plt.title("Average classified pixel value");

Assessment with standard deviation¶

We can also see if the standard deviation of each pixel in time is a reasonable way to learn about the landscape. Similar to how we calculated and plotted the mean above, you can calculate and plot the standard deviation by using the std function in place of the mean function. The standard deviation gives us an idea of how variable a pixel has been over the entire period of time that we looked at.

If you’d like to see the results using a different colour-scheme, you can also try substituting cmap="Greys" or cmap="Blues" in place of cmap="viridis".

In [31]:
plt.figure(figsize=(8, 7))
ds.water.std(dim="time").plot(cmap="viridis")
plt.title("Standard deviation of classified pixel values");

Practice now¶

Make a copy of this notebook so that you can keep this one as a reference. Some things you might like to try:

  1. Find a different area containing water in Africa and re-run the analysis.
  2. Try out a different automatic threshold. How much does it impact the results? See scikit-image for more details about different thresholding methods